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We review a novel approach to evaluate the ground-state properties of many-body lattice 
systems based on an exact probabilistic representation of the dynamics and its long 
time approximation via a central limit theorem. The choice of the asymptotic density 
probability used in the calculation is discussed in detail. 
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1. Introduction 

The real- or imaginary-time dynamics of systems described by a finite Hamilto- 
nian matrix, representing bosonic or fcrmionic degrees of freedom, admits an exact 
probabilistic representation in terms of a proper collection of independent Poisson 
processes 1 > 2 > 3 . For a lattice system, the Poisson processes are associated to the links 
of the lattice and the probabilistic representation leads to an optimal algorithm 
which coincides with the Green Function Quantum Monte Carlo method in the limit 
when the latter becomes exact 4 . 

In the recent Ref. 5 we have exploited the above probabilistic representation to 
derive analytical expressions for the matrix elements of the evolution operator in 
the long time limit. In this way, the ground-state energy as well as the expectation 
of a generic operator in the ground state of a lattice system without sign problem 
are obtained as the solution of a simple scalar equation. The result is based on the 
application of a central limit theorem to the rescaled multiplicities of the values 
assumed by the potential and hopping energies in the configurations dynamically 
visited by the system. As a consequence, the probabilistic expectations can be cal- 
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culatcd by using a Gaussian-like probability density. In this paper, we briefly review 
the approach developed in Ref. 5 and discuss in detail the choice of the asymptotic 
probability density used in the calculation. 

2. Exact probabilistic representation of lattice dynamics 

We illustrate our approach in the case of imaginary-time dynamics for a system of 
hard-core bosons described by the Hamiltonian 

k = - E E%4^+v, a) 

where A C Z d is a finite d-dimensional lattice with |A| sites and Ci a the commuting 
destruction operators at site i and spin index a with the property cf CT = 0. The po- 
tential operator V is arbitrary, e.g. for the Hubbard model V = SieA^ iC iT c iT c ll c ii' 
For simplicity, we assume r\\j — e if i and j are first neighbors and rjij = otherwise. 

In order to study the ground-state properties of the Hamiltonian H it is 
sufficient to evaluate the long time behavior of J2 n ( n \ e ~ Ht \ n o)' where n = 
(nif, nil, ■ ■ • j n |A|T' n |A||) are the lattice occupation numbers taking the values 
or 1. In fact, the ground-state energy is given by 

E = lim -d t logV(n| e -^|n ), (2) 

t— >oo * — J 

71 

while the quantum expectation of a generic operator O in the ground state of H 
can be obtained via the Hcllman-Feynman theorem 5 by evaluating the ground-state 
energy Eq(£) of the modified Hamiltonian H + £0. 

At any finite time t, the matrix elements of the evolution operator considered 
above admit the exact probabilistic representation 

^(n|e-^K) = E(M^ ), (3) 

n 

where is a stochastic functional defined in terms of independent Poisson pro- 

cesses associated to the links of the lattice, see Ref. 3 for a detailed description. At 
each jump of a Poisson process relating sites i and j with spin a and taking place 
at a given configuration n, a particle of spin a moves from site i to site j or vice 
versa if the mod 2 sum of the occupations of these two sites is Xij a (n) = 1, while 
the lattice configuration n remains unchanged if Ay (T (n) = 0. Hereafter, links with 
\ja — 1 w iU be called active. By ordering the jumps according to the times Sk, 
k = 1, . . . , Nt, at which they take place in the interval [0, t), we define a trajectory 
as the Markov chain n\,ri2, ■■■ ,n^ t generated from the initial configuration tiq. 
The number of jumps N t is, of course, a random integer associated to each tra- 
jectory. We associate to each trajectory also two sequences, Aq, A\, . . . , Aw t _i and 
Vo, Vi ■ ■ . , Vfj t , representing the number of active links and the potential energy of 
the visited configurations 

Ak = E A y>(™fc)> ( 4 ) 
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V k = (n k \H\n k ) = (n k \V\n k ). (5) 

Here, T is the set of system links, i.e. the pairs with i < j and i,j G A 

such that rjij ^ 0. The stochastic functional M} nQ which appears in Eq. (3) actu- 
ally depends on the jump times s\, S2, ■ . ■ , s./v t and on the corresponding sequences 
A ,Ax,..., Aiv t _i and V ,Vi..., V Nt . 

3. Probabilistic expectation in the long time limit 

Evaluating the expectation E (.M* ) over the detailed sequences above specified can 
be done numerically by a Monte Carlo method 3 . In Ref. 5 we have demonstrated 
that an analytical expression of E (.M* ) can be obtained in the limit of long times. 
This result is reached in four steps described in the next subsections. The crucial 
point is that, if one integrates over all the possible jumps times, what matter are not 
the detailed sequences Aq, A\, . . . , Ajq t -i and Vq, Vl . . . , Vjf t but the multiplicities 
Na and Ny of the possibles values which the variables A and V may assume. We 
call si and "V the sets of these values and and their cardinalities. It is 
clear that the nature of these sets depends only on the structure of the system 
Hamiltonian, not on the values of the Hamiltonian parameters. The expectation 
E (M f no ) is reduced to an average over Na and Ny- For t — ► oo, this average can be 
evaluated analytically by using saddle-point techniques and a central limit theorem. 

3.1. Canonical decomposition of the expectation 

Referring to Ref. 5 for the details, we decompose the expectation as a series of 
conditional expectations with a fixed number of jumps (canonical averages) 

oo 

E (^n ) = Y, E ( M no\ N * = N )- ( 6 ) 

N=0 

Integrating over the N jumps times, each term of the series (6) can be written as 
E (M^IN = N) = /w N (t) l[A N A , (7) 

\ A£a/ I TV 

where {■) N means average over the trajectories with N jumps generated by extract- 
ing with uniform probability one of the active links available at the configurations 
rao,ni, . . . ,njv-i, and Wjv(i), named weight, is defined as 

W N (t) = e N f dst f ds 2 ... f dsNe -Vosi-v 1 (s 2 - Sl )-...-v N (t~ SN )_ (8) 

3.2. Evaluation of the weights 

According to their definition, the weights satisfy a recursive differential equation 
which is easily solved in terms of the Laplace transform Wn(z) 5 , i.e. 

W*(*) = e"II UTVW (9) 
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While this expression shows that WV(i) depends on the multiplicities Ny for any 
value of N , the explicit inversion of the Laplace transform can be done analytically 
only for N large. However, this is the limit we are interested in since the weights 
Wjv(t) have a maximum at some N which increases by increasing t. By using a 
complex saddle-point method which is asymptotically exact for N — > oo, we get 5 

gZot-E V€ -r N v Iog[(!Bo+V)/e] 

VM*) = 1 e2jVv ^' (10) 

y^Ylve'f {x Q +vy 

where Xq is the solution of the equation 



^ Xq + V 



t. (11) 



3.3. Canonical averages via a central limit theorem 

To evaluate the canonical averages it is useful to introduce the frequencies, v v = 
Ny/N and v A = Na/N, which for N large become continuously distributed in the 
range [0, 1] with the constraints 

E v v = E V A = L ( 12 ) 

Note that for TV large we will not distinguish the different normalizations, N + 1 
and N, of Ny and Na, respectively. Equation (7) can be then rewritten as 

x t+N(v,u) 

duV N {v)—= ^ , (13) 

where v and u are vectors with m = my + m^ components defined as v T = 
(. . . Vy v A . . .) and u T = (. . . — log[(xo + V)/e] . log A . . .) , respectively. 

For later use, we also define v T = (. . . (xq + V)^ 1 .. .) and w T = (. . . (xq + 

V)~ 2 .. .). Note that u, v and w depend on v through xq = xq{v) and 

v = -d Xo u, w = -d Xo v. 

The probability density Vn^) is given by the fraction of trajectories branching 
from the initial configuration no and having after N jumps multiplicities Ny = 
VyN and Na = v A N. For TV large, it can be approximated in the following way. 
We rewrite the multiplicities as Ny = J2k=oXv( n k) and Na = J2k=i^A( n k-i), 
where Xv( n ) = 1 if V( n ) = V and Xv( n ) = otherwise, and similarly for \a- 
Since the configurations to& form a Markov chain with finite state space, a central 
limit theorem applies to each rescaled sum Ny /y/~N or Na/Vn 6 . Therefore, these 
rescaled variables are completely described in terms of the mean values V\N and 
of the covariance matrix S, which are easily measured by sampling over trajectories 
with a large number of jumps. Due to the constraints (12), it is easy to see that 

E v v = E v a = !> (14) 



Wiv(t) n ana 



A^ai 
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Y z ay = y ^v, a = Y = = °> °^uf. (15) 

ver v<£-v A6rf A^si 

In the next two subsections we will describe two different choices for the density 
Vsiy). The two densities differ only for the fact that the first one satisfies the 
constraints (12) in mean, while the other one identically. We will show that, at least 
in the limit of large N, the results obtained by assuming the two densities coincide. 



3.3.1. Purely Gaussian density 

In this case, we assume a purely Gaussian probability density 



*~M- f f ^ < E "'"-"-">. (">) 

Due to Eq. (15), we have that dctS = 0, i.e. £ is singular and the density (16) 
ill defined. Nevertheless, this singularity is eliminable and the calculation with this 
density has the advantage to be rather simple. We rewrite Eq. (13) as 

^Wjv(t) n aNa \ = J due N ^R{u), (17) 



where 



<j>{u) = x ^ + (u, u) - \ (IT V - F), (u - vj) (18) 

and R(u) is a smooth function. The parametric dependence of <f> and Ront and TV 
is omitted for brevity. We will perform the integral (17) with a saddle-point method 
by solving the m-dimcnsional stationary problem for <f>(f). The derivatives of <fi with 
respect to the V and A components of v read 

dv v <t>(v) = jf d u v x - (",v)d„ v x + u v - (E-^i/- V)) v 

d^) = u A -(^- 1 (u-v)) A . 

On the other hand, by using the definition of xo(v), Eq. (11), it is easy to see that 
(i>, v) = t/N, so that <j>{v) is stationary for v = v sp solution of the equation 

i/ sp =77+Ew(i/ sp ). (19) 

Note that, due to Eqs. (14) and (15), we have J2v& r v v ~ ^Ae^s? v& a ~ 1 an< ^ on ^y 
m — 2 equations in (19) arc independent. The result of the saddle-point integration 
then is 



W N (t) TT A Na ) =e N ^R(^)J - ,, (2 ! 7 ) " rr , (20) 
/ v J \ N rn det V 2 <b(v s v)\ 



N 



where 



(21) 
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and V 2 4>(i' sp ) = — S _1 — A(v sp ), A being the m x m matrix with elements 

A a .p = ^-, ^efu/ (22) 

By writing det V 2 0(i^ sp ) = - dct S" 1 dct(l + HA) and observing that detS -1 
cancels out, Eq. (26) becomes 

( W N (t) TT A Na ) = 1 / . (23) 

\ 11 l s V|dct(l + SA)| y/2«Ne>{ V ,w) _ sp 

Note that for V = the matrix A is uniform and det(l + SA) = 1. In general, 
det(SA) = so that det(l + HA) ~ 1 up to terms of second order in A. 

3.3.2. Density with constraints identically satisfied 

In the previous subsection we have evaluated the canonical averages by exploiting 
a central limit theorem for the rescaled variables Ny/^J~N and Na/VN- However, 
due to the constraints (12), the joint probability for these m rescaled sums is not 
purely Gaussian. Given an arbitrary set of m~v — 1 V-likc components and — 1 
A-like components, we can assume a joint probability density as the product of a 
Gaussian density for this set of m — 2 variables and two delta functions which take 
into account the constraints (12), i.e. 



(24) 



where Ts(y) is the normal density defined in terms of the vector v having the m— 2 
chosen components of u 



^= f m ^tr [ e-«*-<>-*>.<>-*». (25 ) 

In the following, with the symbol * over a vector or a matrix we will indicate the 
projection onto the chosen m — 2 dimensional index space. 
By using the density (24), we rewrite Eq. (13) as 



n 



A Na ) = / dve N ^R{£>), (26) 



where R is a smooth function of i> and, if V* and A* are the two components chosen 
to reduce to m — 2 the dimension of the index space, 

=*0T7 + (»>.«(£()))+ [1- P V ) u vA £ a)+ [1- X! ^ "A, 



if±-^-F),(i>-P) 



2 



l 



I 
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Note that £q = xq(0) is the solution of Eq. (11) where the constraints (12) have 
been explicited. The function cj) can be rewritten as 



t 



*)=*0T7+ [OM )-- S \v-v),(v-v) 



N 

where Su = u—Uq and Uq = (u v . . . u v ; u A . . . u 
of <j> with respect to the components of v , we get 



(28) 



) . By evaluating the derivatives 



dp. 



) = W^v^o- [i>,5v)dp v x +Su v - (S l (v-v 



d i >k£>) = 8u A -{±- 1 {i>-i>) 



v V, d i> v Xo 



where Sv = —d$ Su = — (d% u — d$ Uo). Observing that (v, Sv) + v v = t/N, we 
can rewrite these derivatives as 

dejiv) = 6u v - (£- l (u-v) 
doJ{0) = 5u A -(t- 1 {v-V) 

Therefore, the saddle-point equation is 



O sp ^ 17 + T,5u(u sp ) 



(29) 



Finally, we calculate <^>(z> sp ). By using Eqs. (28) and (29) and the identities ( v>, Su 



i = (u,u) and ("SSu,Su) = (Sit, u), we get 



(*"P) 



•?'o 



— + (v + IlSu, Su) - - (t,Su, Su 



N 



2 V 



+ u Vm + u At 



(30) 



Due to the identity (^ESu, Sv) = (Sit, v), we have xo(C> sp ) = xq(u sp ), where xo(u sp ) 
is the value obtained with the purely Gaussian density. Therefore, Eq. (30) shows 
that, at least at the saddle-point level, the chosen probability density provides the 
same results obtained in the previous subsection. 



3.4. Resumming the canonical series 

In order to evaluate the expectation E (M f no ) we need to resum the series (6). For 
t — > oo, the sum can be substituted with the integral 



dN 



^|det(l + SA)| J v / 2TrNe 2 (v sp ,w(xo(v sp )))' 



where 

i/>(N) = x (v sp )t + N 



(v,u(xo(» sp ))) + ^(Su(x (^ p )),«(xo(^ p ))) 



(31) 



(32) 
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The integrand in (31) is exponentially peaked at TV = iV sp , where iV sp satisfies 
dNip(N sp ) = 0. For a generic N we have 



{t-N [(17, v) + (Et», «)]} W &„ + (17, «) + - (Su, «) 



where, due to Eqs. (11) and (19), we noticed that the term {t — N \{y, v) + (Six, v)]} 
vanishes for v = v sp . The stationarity condition for ip(N) is, therefore, 



(i/, u) + — (Sit, u) 



0. 



(33) 



i/=v s v ,N=N s p 



Equation (33) is a time independent equation which determines Xo\ v — u ap tv=jv b p as 
a function of V and X. According to Eq. (11), this means that the quantity N sp 
increases linearly with t so that Xq\ v=vSV 7v=jv b p becomes independent of time in the 
limit t — > oo. By evaluating the second derivative of ip(N) 

d 2 N ^(N) = -(is» p ,v(x (v sp )))d N x (» sp ), 



where 



d N x (v sp ) = 



(v sp ,v(x (is sp ))) 



N (u sp ,w(x (v sp ))) + (v^o^p)),^-^^ 815 )))' 



(34) 
(35) 



and approximating ip(N) ~ ip(N sp ) + ^d%ip{N Bp )(N - N sp ) 2 in the exponent and 
N ~ iV sp elsewhere, the integral (31) reduces to a Gaussian one which gives 

/ l + tr(5L4) e Xot 



|det(l + SA)| e(u,v) 



(36) 



Note that, since the peak of the Gaussian at N = N sp moves to infinity linearly 
with t while its width increases only as y/i, the result (36) is asymptotically exact 
for t — » oo. 

According to Eq. (2), the result (36) shows that the ground-state energy of the 
hard-core boson system is 

Eob = - x o\ v= ^ P N=N s P ■ (37) 

Equation (33) is, therefore, the equation for the ground-state energy It defines Eqb 
in terms of V and £ and explicitly reads 



= - £ XVlogj ^L±Z) + £ VAlog(A) 



verve/ 



log 



-E 0B + V 



£ E S ^lo g r^ B + V peg (A) 
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+ 1 E E ^ log (A) log (A'). (38) 
In the case V = 0, the ground-state energy £^ can be solved analytically 



c(0) 

E 0B = - e ex P 



4. Conclusions 



J2v A log(A) + lYl E S^log(A)log(A') 



2 



(39) 



By using saddle-point techniques and a central limit theorem, we have exploited an 
exact probabilistic representation of the quantum dynamics in a lattice to derive 
analytical approximations for the matrix elements of the evolution operator of a 
system of hard-core bosons in the limit of long times. The approach yields a simple 
scalar equation for the ground-state energy. This equation depends on the values of 
the generalized potentials V and of the kinetic quantities A, and on the statistical 
moments V and X of their asymptotic multiplicities Ny and Na- In turn, these 
moments depend only on the structure of the system Hamiltonian, not on the values 
of the Hamiltonian parameters. This implies that the statistical moments must be 
measured una tantum for a given Hamiltonian structure and, once V and S arc 
known, our approach provides the ground-state energy analytically as a function of 
the Hamiltonian parameters. 

The ground state energies obtained with the present formula reveal a small 
systematic error when compared with exact results 5 . In fact, a more accurate 
analysis which takes into account the large deviations neglected by the central limit 
theorem shows that the present approach corresponds to a second order truncation 
of an exact cumulant expansion for the ground-state energy 7 . 

Similar results hold in the case of Hamiltonians with arbitrary kinetic operators. 
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